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Abstract: In emission tomography, the process of image reconstruction has 
evolved since the utilization of the filtered back projection algorithm (FBP). 
Nowadays, iterative techniques incorporating statistical knowledge of the pho¬ 
ton emission process has produced much interest among researchers, not only 
for the improvements achieved in image quality, but also by the possibility of 
incorporate corrections factors, like scattering, attenuation, etc. 

One of the problems that has not been discussed in depth is the incorpo¬ 
ration of correction methodologies due to patient breathing, which produces 
blur regions in lungs and cardiac images. Solutions like respiratory gating, 
that synchronize the breathing cycle of the patient and the data acquisition 
process, or correlated dynamic PET techniques that use external radioactive 
markers and list mode data, have been tested with improvements over the 
spatial activity distribution in lungs lesions, but with the disadvantages of 
requiring extra hardware or more expensive scanner systems. The objective 
of this study was to incorporate breathing-movement corrections directly to 
the phase of image reconstruction, without any additional acquisition protocol 
consideration. For this, a procedure of correction inside the probability ma¬ 
trix in the classical MLEM algorithm has been implemented which takes into 
account the spatio-temporal relationship of each voxel in the structure under 
study. We present 2D and 3D results from synthetic and realistic simulations 
showing the potential benefits of using this approach. 
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Correction du mouvement respiratoire en 
tomographie d’emission 

Resume : En tomographie d’emission, les precedes de reconstruction ont 
evolue des methodes analytiques de retro-projection filtree, aussi utilisees en 
tomographie de transmission, a des techniques iteratives qui incorporent une 
modelisation statistique de remission des photons par un element radioactif. 
Ces dernieres ont suscite un grand interet dans la communaute, non seulement 
par leur aptitude a ameliorer la qualite des images reconstruites, mais aussi 
par leur capacite a incorporer divers facteurs de correction, prenant en compte 
la dispersion, l’attenuation, etc. 

Toutefois, la prise en compte des mouvements du patient, en particulier des 
mouvements physiologiques (respiration, battement cardiaque), qui introduis- 
ent un flou des l’acquisition des projections, flou que Ton retrouve lors de la 
reconstruction, n’a ete aborde que superficiellement. 

Jusqu’a present, les techniques employees visent a adapter l’acquisition des 
donnees pour rendre la reconstruction d’un objet statique possible. II s’agit par 
exemple de synchroniser l’acquisition des donnees avec le mouvement (comme 
la respiration), ou d’acquerir les detections en mode liste et d’utiliser alors 
une technique de correlation dynamique qui necessite l’ajout de marqueurs 
radioactifs externes supposes rendre compte du mouvement. 

Dans le cadre de l’imagerie des poumons, ces techniques ont montre une ameli¬ 
oration de la quantification des lesions detectees, mais malheureusement au 
prix d’un protocole d’acquisition plus lourd (necessitant des materiels supple- 
mentaires et/ou plus couteux), qui les rendent inapplicables a posteriori a la 
plupart des acquisitions en routine clinique. 

L’objectif de ce travail est d’incorporer directement les corrections lies au 
mouvement dans le processus de reconstruction, sans avoir a agir en amont 
sur le processus d’acquisition. Pour ce faire, la technique proposee s’appuie 
sur un algorithme MLEM (Maximum Likelihood, Expectation Maximization ) 
dans lequel la matrice des probability prend en compte la contribution spatio- 
temporelle a chaque projection du voxel a reconstruire. 

Nous presentons des resultats 2-D et 3-D, avec des experiences synthetiques 
et des simulations realistes, qui montrent le benefice potentiel d’une telle ap- 
proche. 
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1 Introduction 

Nowadays, several medical imaging modalities exist that allow the physicians 
to get an image-based description of an organ. This description can be ei¬ 
ther information of the anatomical structures under study or information of 
their functional behavior. This functional information has brought interest 
among physicians, physiologists, biologists, etc. due mainly to the richness 
of the information in the sense of detection and treatment of diseases like 
cancer in which the biochemical changes precede the anatomical ones, allow¬ 
ing earlier diagnosis and treatment. Furthermore, the combination of func¬ 
tional and anatomical modalities has also allowed improvements in diagnosis 
and following of diseases. All this, in addition with advances in data ac¬ 
quisition (i.e., higher resolutions, lower acquisition times, etc.), has caused a 
demand of better algorithms, like registration techniques to fusion modalities 
[40, 30, 29, 28, 47, 48, 7, 46, 14, 8, 9, 25], or image reconstruction algorithms 
that take a set of indirect measurement data and generate a visual representa¬ 
tion of a variable of interest. For the latter, the widespread improvements in 
the techniques include better image resolutions, inclusion of physical correc¬ 
tions, faster reconstruction algorithms, etc. [18, 26, 1, 43, 36, 27] 

Unlikely to CT, MRI and others techniques that use an external source of ra¬ 
diation to visualize the different structures (measuring absorption coefficients 
(CT), protons density (MRI), etc.), in emission tomography (EM) an inter¬ 
nal source of radiation (commonly injected) is used, called radio-isotope. The 
type of radio-isotope makes a difference between EM techniques, being PET 
(Positron Emission Tomography) and SPECT (Single Photon Emitting Com¬ 
puted Tomography) the classical ones. The radio-isotope goes "attached" to a 
tracer which is selected depending on the target organ to be visualized. This 
process is not free of artifacts, the principal factors are: body attenuation, 
scattering, hardware factors, etc. Several studies have shown different ways to 
correct for body attenuation, scattering and so on, but little effort has been 
done for the case of patient movement correction. For example, in the case of 
lungs studies, the breathing cycle affects the quantification of lesions. In [32] 
the authors investigated the impact of movement during the exam, finding 
mislocalizations of lesions in the fusion of PET and CT. In [34] they found 
up to an 34% of lesion volume reduction. Significant tumor motion have been 
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reported in similar studies [44, 45, 3]. Some studies have been performed in 
order to correct for breathing movement, they include respiratory gating [34], 
which synchronize the breathing cycle with the data acquisition process in or¬ 
der to get only data of a specific phase within the breathing cycle, others use 
respiratory-correlated dynamic PET, where an external FDG source point sit¬ 
uated over the patient’s thorax correlates with the patient’s breathing motion 
[35]. These techniques have shown improvements in reducing the smearing ef¬ 
fects caused by the breathing movement, contributing to a better quantification 
of lung lesions, in the other hand they either require extra hardware (respira¬ 
tory gating) or specific data acquisition modes like list-mode data, where not 
always is possible to have access to this type of data acquisition mode. In 
this article we provide a methodology that works directly in the accumulated 
projection set (sinogram), correcting for breathing movement. The method 
have been tested with phantom images produced with the SimSET simulator 
and the thorax phantom NCAT. 
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2 Image Reconstruction Algorithms in Emission 
Tomography 

As it was said, the aim of an image reconstruction algorithm is to obtain from 
a set of indirect measurements a visual representation of a physical quantity. 
In Emission Tomography (EM), this physical quantity corresponds to the spa¬ 
tial radiation density, and the indirect measurements correspond to photon 
counts for every angle, recorded by scintillation detectors. Once the counts 
are acquired, (i.e., projections) an algorithm is applied to obtain the spatial 
distribution of the radiation density generated for the radio-isotope which was 
injected into the patient’s bloodstream or inhalated. 

These algorithms can be classified in two types; analytical and algebraical. We 
present a brief description of the most common algorithms that is by no means 
a complete description of them, but just intended to show a general view. 

2.1 Analytical Algorithms 

The analytical algorithms are based in an analytical model of the process of 
acquisition, this model is based on the Radon transform 

/ OO /»oo 

/ f(x, y)5(xcos0 + ysinO — t)dxdy (1) 

-OO J — OO 

Where m(t , 6) represents the set of line integrals passing through the object 
f(x,y). Thus, the problem consists in find m(t,9 ) from f(x,y). Application 
of inverse FFT doesn’t work well because it implies interpolations at high fre¬ 
quencies where the density of the 2D Fourier space is low [1]. The most typical 
inverse Radon transform algorithm is the Filtered Backprojection algorithm 
(FBP). It consists of taking each projection profile, take its Fourier Transform 
(FT), apply a frequency domain filter, take the inverse FT and backproject 
over the image at the given angle. In emission tomography due to the noisy 
nature of the data and the incompleteness of the projection data, images re¬ 
constructed with FBP suffer from heavy noise. Furthermore, the filtering step 
amplifies high frequencies, incrementing the noise level. This effect can be 
decreased applying a window function to the filter but with the need of con¬ 
sidering a trade-off between noise level and image resolution due to overlapping 
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of frequency bands between the true and noise signal. FBP is fast and not hard 
to implement. However, it oversimplify the data acquisition process consider¬ 
ing all measurements equally, which leads to noisy images. In the next section 
a review of iterative algorithms is presented. 

2.2 Iterative Algorithms 

Contrary to analytical algorithms, iterative algorithms allow better modelliza- 
tion of the acquisition and emission process. Also, the modelling is discrete 
and not continuous as in the analytical case (not considering the impleme- 
mentation). Generally speaking the idea behind consists in, given a set of 
measurements p and the projection matrix R, find the set of values / that 
accomplish the relation p = Rf. It will be shown later how the matrix R can 
be constructed and what other type of information can be added to it. 

Use of direct algebraical methods to obtain / is not possible due to the large 
size of the matrix R. Besides, noise in p and the approximation of R doesn’t 
allow an exact solution of / [1]. Furthermore, use of least-squares and pseudo¬ 
inverse may deliver negative values. 

Algebraical methods overcome these problems using an iterative approach. For 
every iteration, a projection of the guess image is performed, which is com¬ 
pared, by means of some criteria, with the measured data (i.e., p). The error 
produced is feedback into the guess image, and a new iteration is performed. 
The problem with this type of approach is that as the iterations increase the 
image becomes noisy. This is because the algorithm searches to get / as close 
as possible to p , which contains noise data. In fact, what is of interest is not 
the complete data / (true data + noise) but its mean A. In the next section 
we review the Maximum-Likelihood Expectation Maximization (MLEM) algo¬ 
rithm which tries to find the mean values of / by means of modelling the data 
acquisition as a Poisson process. 
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2.3 Maximum Likelihood Expectation Maximization Al¬ 
gorithm 

This algorithm, first introduced in emission tomography by Shepp and Vardi 
[43], is based in a Poisson model for the emission process. For a given emission 
element b the number of emissions f b follows a Poisson law with mean \ b 

f b ~ Poisson(Xb) 

Besides, the projection matrix R (or called by some authors system matrix or 
transition matrix ) gives the probability that a certain emission from pixel b is 
detected by the detector d (called dexel hereafter). Furthermore, the number 
of detections from dexel d (i.e., pd ) can be calculated as the sum of detections 
coming from every emission element b 


b=n 

Pd = ^Pdb ( 2 ) 

b= i 

Furthermore, the number of detections from emission element b detected by 
dexel d can be written in function of the number of emissions f b 

Pdb = fbRdb (3) 

From (2) and (3), the number of detections can be expressed in terms of the 
number of emissions, as follows 


b=n 

Pd = Yl f bRdb ( 4 ) 

6=1 

Similarly, for the mean counts we can write 


b—n b=n 

Pd -^[Pd] ^ ^ Xdb ^ ^ X b Rdb (5) 

6=1 6=1 

Where Xdb stands for the mean number of emissions from emission element b 
being detected by dexel d. 
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The ML estimate searches iteratively to find a value A that accomplish the 
following 

A = arg max [/(A)] (6) 

A 


with /(A) 


/(A) = P(p |A) 


( 7 ) 


the likelihood of getting a set of measures p given the image A. 


To simplify calculations, it is common to employ the log-likelihood, L(A) = 
log[/(A)]. More important is the fact that the variables A db are independent 
variables. Thus, L(A) can be expressed as: 


L{ A) = log(P(p/A)) = log 



e -A d b^l\ 
Pdb \ ) 


And equation (6) becomes 


A 


arg max [L( A)] 

A 


( 8 ) 

(9) 


In order to solve for (9), maximization of (8) is performed. Taking first deriva¬ 
tive of 8 


dL(X) 
d\ b 




- - ^ R db+Y^ 

d d 


Pdb 

A b 


0 


9LW = n t = ZdPdb 

d\b b Ed R db 

However, since we don’t have direct access to pd b , it is necessary to calculate its 
conditional expectation given the measures pd and the current guess image A. 
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Further, remembering that for independent Poisson variables X,Y with means 
Ax, Ay, the expectation of X conditioned on the sum A" + Y is E[X\X + Y] = 
^Ax+Ay Y [43|- The expectation of p db can be calculated as 


E[p db \p d ,X] = 

X,b' A db' 


PdX b R d b 
'Yhb' ^b'Rdb' 


( 11 ) 


Then, substituting equation (11) into (10) we get 


\ <iC+l> _ 
A b ~ 


A^ A > PdRdb 


( 12 ) 


Where the index K stands for the iteration number. 

From (12) we can observe that the algorithm basically follows the same prin¬ 
ciple already mentioned for the iterative algorithms. It takes a guess image 
\ <K> , makes a re-projection of it (denominator in right summation). Then, 
the ratios between the measured data and the forward-projected values are 
backprojected into the image space and added, the sum of this values are 
used as multiplicative update coefficients to iteration K + 1. The term J2 d R d b 
acts as a normalization value, as in theory J2 d R d b = 1, meaning that every 
emission should be detected by the system. 


The principal characteristics of the MLEM algorithm are its non-negativity 
(i.e., it assures non-negative pixel values for all the images generated) and, for 
every iteration the number of emissions equals the number of detections. 

One aspect that still remains open is when the iterations should be stopped. 
Several measures exist that can be used to check the quality of the recon¬ 
structed image. 

It has been shown [26, 21] that a good figure-of-merit is the Root Mean Square 
value (RMS) 


RMS k = (13) 

V h 

The RMS value shows that in a point in the iterations, its value reaches a 
minimum and then starts to increase, indicating that noise in the measured 
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data begins to be added to the reconstructed image [22], 

While the RMS value has only utility in simulations studies, where the density 
distribution can be known a priori, it can not be applied to real studies. 

The likelihood of the objective function can be used as a statistical criteria to 
stop the algorithm: 

L W = Y^ Pd log {pa) ~Yd~ log(p d !)] ( 14 ) 

d 

Using equation (5) in (14), the likelihood can be calculated as: 

L( A) = Y [Pd log ( y A b R d b | - Y ^Rdb ~ log(p rf !)] (15) 

d \ b J b 

The problem of using equation (15) is that as the iterations continue, the like¬ 
lihood will increase (monotonicity of the solution), without indicate the point 
where noise will begin to be added to the reconstructed image. 

In [10], the author presents a new stopping criteria, that consists in sepa¬ 
rate randomly the projection data in two halves, namely A and B. Then, one 
proceeds with the reconstruction using the set A, and for each iteration the 
likelihood of the data set B is calculated using the estimates obtained with 
A. It is shown that the likelihood will increase until a certain point where the 
iterations over A are stopped, changing to the data set B. Once both points 
of convergence are reached, the two estimates are summed up to obtain the 
final image estimate of the density distribution. This technique has shown 
good results on noise rejection but as it has been remarked in [19], the cross¬ 
likelihood is dependent of the number of counts. In [19], Johnson proposed 
a variant to the cross-likelihood scheme, in which the projection data set is 
divided in k subsets. Then, each subset is subtracted from the complete data 
set. Each subtracted data set is reconstructed and multiplied by 1 /(k — 1) 
to preserve the counts number in each iteration. For each subtracted data 
set, the iterations are stopped when the likelihood of the non-included data 
set is maximized. Other approach has been proposed in [22], where a study 
of the multiplicative update coefficients of the MLEM algorithm has allowed 
the authors to establish a stopping rule. For each iteration the update coeffi¬ 
cients are stored and histogrammed. The technique is based in the fact that 
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the optimum iteration value (given in a simulation study by the RMS value) 
is reached always in the same value for the histogrammed coefficients. In a 
precedent work [23], the authors stated that a value of 0.8 produces images 
close to the optimal reconstructed image (with ±5 iterations). 

The MLEM algorithm has shown better results than the classical FBP algo¬ 
rithm, however its inconvenient is its slow convergence. Different approaches 
has been created to overcome this problem, like parallelism of the process, com¬ 
putational memory management improvements to speed up iterations, and so 
on. One technique that has improved the speed of convergence and has ap¬ 
proached the EM technique to the clinical scenario is the Ordered Subsets 
Expectation Maximization algorithm (OSEM) [18]. Instead of working with 
one set of projections, the algorithm performs several sub-iterations over a 
smaller subset before beginning with the next one. The results using this tech¬ 
nique has shown that the image reconstructed gets closer to the convergence 
point in less time than the original MLEM algorithm. It occurs because for 
each iteration, a voxel is updated as many times as number of subsets. There¬ 
fore, each voxel will be visited more times for iteration than in the case of 
using a single set of projections. In other terms, for every subset an image is 
reconstructed considering the information contained in that subset, then, for 
the successive subsets new information is added to the reconstructed image. 
The speed up of the algorithm is based on the smaller time required to recon¬ 
struct an image with less projection data. According to [18], the higher the 
number of subdivisions the better level of detail it can be obtained. However, 
they agree that there is a limit with the number of subdivisions. Beyond this 
limit, the algorithm lacks of sufficient data to fit the observed data. 

Other issue of interest is how to select and how to order the subsets. In [18], 
the authors remark that the selection of subsets should be done in a balanced 
way, so that the voxel activity information is also balanced in the subsets [18]. 
Regarding the order in which the subsets are processed, They suggest that 
even if the order is arbitrary it is preferable to process subsets which include 
new information as soon as possible. 

Despite the good empirical results presented by OSEM, its major problem 
has been its lack of proof convergence and ambiguity in the use of priors [1], 
also it has been stated that it can lead to limit cycles in the iterative object 
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estinates [17]. Recently, an accelerated convergent ordered subset algorithm 
was presented [17], that establishes a tradeoff between speed and convergence 
by using a parameter that updates itself automatically as the iterations pro¬ 
ceeds. This parameter introduces a linear combination between the fast but 
non-convergent OSEM algorithm and the slow but convergent COSEM algo¬ 
rithm [27]. Basically, at the beggining of the iterations more weighting is given 
to the OSEM image, in order to speed-up the convergence, and then, as the 
iterations proceeds, more weighting is given to the COSEM image in order to 
ensure convergence of the global algorithm. The major drawback is that the 
linear combination needs a precalculation of each guess image (i.e. OSEM and 
COSEM), which slows the algorithm. The authors stated that this is a point 
of further research. 


3 Incorporating movement correction into the 
image reconstruction algorithm 

How it was said, the correction of movement has been performed mainly by 
using components attached to the data acquisition system, in an effort of 
neglecting information that disagree spatially with a static reference, like a 
CT image, where a fusion of modalities is wanted. Differently, our approach 
is based in the information of the breathing cycle that can be included in the 
step of image reconstruction. 

Different authors already noticed the convenience of introducing correction 
factors into the system matrix, like attenuation factors, scattering factors, 
detector efficiency, etc. This convenience comes from the fact that each element 
of the system matrix corresponds to the relationship between a given emission 
element and a detector element. This one-to-one relationship permits the 
inclusion of these factors in a straightforward way. 

To correct for respiratory movement, we have introduced a correction factor 
that takes into account how an element of emission moves and deforms itself 
within the respiratory cycle. 

In the case of no movement during the data acquisition period, each element 
of emission will contribute with emissions to a given set of dexels. But, in the 
presence of movement, it is more likely that the quantity of detections that 
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each dexel receives will change (even to the point where some dexels will not 
receive any contribution from the given emission element). Considering this, 
and remembering that each element Rdb of the projection matrix represents the 
probability that an emission from b is detected by dexel d. allow us to correct 
for movement by recalculating the values of the projection matrix by taking 
into account the position of every emission element during the data acquisition 
period. Let’s model the movement as a discrete set of time-weighted positions. 
Then, from equation (12) we can define R db as the detection probability of an 
emission produced by the emission element b and detected by dexel d during 
the movement state k. Therefore, the new detection probability is: 

Rl = £ w, (16) 

i 


Where each Wi(i = 1,..., W) is a time weight for the state i and W is the total 
number of discrete positions for our discrete movement model. 

The time-weight takes into account the fact that the more time a moving 
emission element takes over a detection element, the more number of detections 
this element will capture. 

Finally, introducing (16) into (12) we get 


\ <K+1> 
A b 


K K> PdRdb 

E.-RS, £* x i K>Rl z* 


(17) 


Equation (16) tries to relocate the number of emissions that due to movement 
were detected in others dexels but the dexel where they should have been 
detected in no presence of movement. 

The way each term is calculated in equation (16) will take into account the 
position of the emission element at each state of movement and is calculated 
as follows 


R l db = 


l db 


ld'b 


(18) 


Where the term P db is the intersection length between the dexel d and the emis¬ 
sion element b during movement state i. The summation in the denominator 
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Figure 1: Calculation of system matrix elements considering movement, 
of equation (18) acts as a normalization term. 

Fig. 1 schematize the way the movement of an emission element will affect the 
calculation of equation (18). Here, at each grid position, each emission element 
bn has been modelled as either a circle (2D) or a sphere (3D). Then, in the 
case of non-movement, the probability R.d,b n 7^ 0 and Rdjb n = 0. But if move¬ 
ment comes into play, the situation is the opposite (Rdib„ = 0 and R. ( i :i b, n ^ 0). 
Furthermore, as remarked in Fig. 1, not only a displacement of the emission 
element is been considered but also deformation of it. 

3.1 Calculating emission-element volume deformation by 
diagonalization of Jacobian matrix 

As it was said, for the calculation of each system matrix term, not only dis¬ 
placements of each emission element are considered but also volume deforma¬ 
tion of it. 

Let’s be U(x,y, z ) the displacement vector field having the information of how 
an emission element with coordinates (x, y, z ) moved in space to the position 
(x p ,y p ,z p ), and ip(x,y,z) the deformation function defined as: 

x p = x + U x (x, y, z) = ip x (x, y, z) 

y p = y + u y {x,y,z) = <Py{x,y,z) \<p(x,y,z) (19) 

z p = z + U z (x,y,z) = tp z {x,y,z) ) 

Then, it is possible to calculate the gradient of p{x } y 1 z) at a certain emission 
element p = (p x , p y , p z ) 


RR n° 5279 



















18 


M. Reyes, G. Malandain, J. Darcourt 


\7 p ip = 


/ 

dip x 

dp x 

d<fx \ 


dx 

dy 

dz 


dpy_ 

dip^ 

dp>y 


dx 

dy 

dz 

V 

d^>z 

dip z 

dip z 

dx 

dy 

dz / 


X = p x 

y = Py 
z = p z 


( 20 ) 


Since V p <p = Id + VU(x,y, z) with Id the identity matrix, equation (20) can 
be written as: 


V p(p = 


/i + f- 

9Uy 

dx 

8Uz 

dx 


9Ux dU x \ 

dy dz ’ 

1 i 9_Uy dJR 

1 “ r dy 
dU 2 
dy 


OUz 1 

a,. 


dz 
I dUz 
-r dz 


X = Px 

y = Py 
z = p z 


( 21 ) 


The value of the determinant of V<^(a.k.a jacobian of p>) describes if the emis¬ 
sion element suffers an expansion (| V p p > 1), a contraction (|V P <^| < 1) or if 
it preserves its volume (|V P </?| = 1) [37]. Moreover, it is possible to calculate 
in what direction and magnitude the emission element will either expand or 
contract. 

Let’s assume that the matrix W p p is a real symmetric matrix. Then, there 
exist matrices P and D, such as: 


V p p = PDP~ l (22) 

Where, D = diag(X^\ with X J p (j = 0,1,2) eigenvalues of W p p and 

each column of matrix P, an eigenvector of matrix V p (p. What is of our inter¬ 
est, is the information that can be retrieved from matrices P and D. 

Since the eigenvectors represent the preferred directions of a matrix, in our 
case they represent the preferred deformation direction of an emission ele¬ 
ment, whilst every eigenvalue gives the magnitude of the deformation in the 
direction of the associated eigenvector. 

Fig. 2 presents a test performed to check how emissions elements (modelled as 
spheres) translate and change their volume according to a pre-built displace¬ 
ment vector field. 
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Figure 2: Testing the deformation of a set of sphere-shape-modelled emissions 
elements following a pre-built displacement vector field. Left: Original set and 
displacement vector field. Right: Translated and deformed set. 

4 Implementation 

This section presents the implementation issues concerning the parallelization 
of the MLEM algorithm as well as the accelerators for the forward and back¬ 
ward operators. To finalize with a brief description of the SimSET library and 
the NCAT thorax phantom, utilized for the lungs simulation. 

4.1 Parallelization of the MLEM algorithm for 3D image 
reconstruction 

In section 2.3, it has been remarked the slow convergence of MLEM, and some 
solutions were presented in order to decrease the convergence time. While in 
two dimensions the speed problem of the technique is feasible to solve by using 
for example smaller projections subsets at each iteration (OSEM approach) 
or by designing more complex ways to store the system matrix to gain faster 
data access times. In three dimensions it is required a force brute method like 
parallelism. 

Parallelization of EM algorithms has been already used to speedup the con- 
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vergence, in [21] a summary of previous works is described. Here the common 
denominator has been the distribution of the computational load of forward 
and backward projections over a set of processors. 

The main aspect of the algorithm in terms of parallelization, is without a doubt 
the independence between neighbor emissions elements over a single iteration, 
which eases its implementation. 

Our implementation is composed by three stages of parallelism. First, as 
suggests equation (16), calculus of the normalization term for the corrected 
detection probability should be yielded before the steps of forward and back 
projection. Then, parallelization of forward projection step is performed. The 
counts estimates are stored and used in the next and final step of paralleliza¬ 
tion; backward projection and image update. 

Fig. 3 shows a diagram where it can be seen how each one of this paralleliza¬ 
tion tasks relates to the others. Each box in Fig. 3 represents an entity of 

sinogram space 
image space 



-Calculated data goes through master to final destination 

_► Direct communication slave/master 

send of slave-specific information 

->- Direct commnicatin master/slave 

send of common information 

Figure 3: Implementation diagram of parallelized MLEM algorithm. 

execution, starting by the master that makes the light calculations (i.e., cal¬ 
culating the sinogram and image space dimensions, reading the sinogram data 
file, initializing the image to be reconstructed, etc.). The first slave created 
is SlaveLS, which calculates the normalization term in equation 16. Once the 
slave has finished its task, it sends the data back to the master, which store 
it to be used later for the two other slaves. The second slave, SlaveFP, is in 
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charge of the forward projection part (i.e., right denominator in equation 17 
). Besides from the basic information delivered by the master, it receives the 
already calculated normalization term. The detection count estimates are sent 
back to the master, which will spawn the last set of slaves that perform the 
back-projection operation. Once again, the data calculated in the previous 
steps is sent to every SlaveBP. Once the SlaveBP has finished its work, it 
sends the portion of the updated image estimate to the master, which gets the 
data from each slave and assembles the image. At this point an iteration is 
completed. 

For the communication between master and slaves, the PVM (Parallel Virtual 
Machine) software was used [11] and execution was performed in the INRIA 
Sophia Antipolis cluster system [15]. 

4.2 Accelerating the MLEM using geometrical consider¬ 
ations 

Geometrical considerations can be used to diminish the number of calculations 
in equation (12). Specifically in the two more computationally expensive op¬ 
erations: forward and backward projection. For the former, the first approach 
was to calculate a region of interest that takes into account the geometrical 
space traversed by a specific dexel (see Fig. 4). Later, the well-known Bresen- 



Figure 4: Region of interest to reduce calculation times in forward projection, 
ham’s algorithm [5] [50] used to represent a continuous line in a discrete space 
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was incorporated as a better accelerator. While the goal of the Bresenham’s al¬ 
gorithm is to better represent a continuous line over a grid space, sometimes it 
fails in adding all the pixels (2D case) traversed by the line. Which in our case 
is of vital importance in order to take into account all voxels contributing to a 
specific dexel. Then, for each voxel given by the Bresenham’s algorithm (3D 
version), a 6-neighborhood was considered to ensure that the voxels traversed 
by the dexel are present. Fig. 5 shows the result from the implementation of 
the 3D Bresenham’s algorithm with the previously mentioned 6-neighborhood. 



Figure 5: Bresenham 3D with 6-neighboring to accelerate forward projection. 
The red cubes, are the ones from the 3D bresenham’s algorithm, and the others 
are the 6-neighbors for each of the red cubes. This, to ensure that all voxels 
traversed by the dexel are included. 

To accelerate the backward projection step, at every voxel an angle-specific 
forward projection operation is performed to get a list of dexels passing trough 
it. Formerly, for a given emission element b , the non-accelerated backward 
projection operation consists on calculate: 



(23) 
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With D the set of all dexels forming the scanner system and n ( j the current 
estimated value for emissions detected by dexel d. As not all dexels d will 
traverse emission element b (yielding null probabilities Pdb). It is possible to 
reduce the set of dexels to be visited by applying a forward projection operation 
to the emission element b at every azimuthal angle O^i = 1, 

D' = [F„. (6),. ..,F „„(!,)] (24) 

Where, F 0 (b) stands for the forward projection operator applied to emission 
element b at angle 0. Finally, the backward projection operation of emission 
element b over the reduced set D' will be 

v-wn , efl , 

V ** 

The previous procedure is depicted in Fig. 6 where continuous lines represent 
dexels in reduced set D' and non-continuous lines are examples of dexels being 
discarded from set D'. 



Figure 6: Reducing dexel space in backward projection. Continuous lines 
represent the dexels that are being considered for the forward projection while 
non-continuous lines are not included in the forward projection step. 

If movement correction will be incorporated, it is necessary to include it in 
the acceleration steps. As it was said for the forward projection step, the Bre- 
senham’s algorithm was used to obtain a speed-up. This speed-up is basically 
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due to the smaller set of emission elements used in the estimation of outcomes. 
In the presence of movement, the emission elements contributing to a specific 
dexel during movement have to be added to the set of emission element given 
by the Bresenham’s algorithm. 

Formerly, for dexel d, the set B of emission elements to be considered during 
forward projection is: 

B = Bresenham(d) u|^J-U(<i),i = 1 , - -,W 

i 

Where Bresenham(d) is the set of emission elements given by the Bresenham’s 
algorithm and 


p\d) = * 0 } 

Similarly, for the back-projection step the presence of movement should be 
taken into account. Revisiting equation (24) and (25), we can recalculate the 
set of dexels passing through a specific emission element by considering the 
states of movement. Then, equation (24) can be rewritten as: 

d' = (JKU'W)..... (&))],; = o,.... w (26) 

i 

with T l {b) the transformation that gives the position of emission element b in 
movement state i. 

The ML-EM algorithm was implemented and tested with synthetic images 
generated with the SimSET library, which is briefly described in the next 
section. 

4.3 Using the SimSET library to simulate emission to¬ 
mography 

The Simulation System for Emission Tomography^ SimSET), is a library that 
simulates the process of emission tomography allowing the configuration of 
several components, going from the physical to the instrumentation process. 
Fig.7 shows these components. 

The principal modules are the Object Editor , in which is possible to create 
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Figure 7: The SimSET library modules. 


activity and attenuation objects (using basic geometrical forms). The Photon 
History Generator , that creates and follows the photons through the different 
layers of the simulation. Finally, The tracked photons can be passed to the 
Binning Module which generates the sinogram data file (among others for¬ 
mats). 

Through a command line it is possible to generate the activity and attenu¬ 
ation objects. In addition, the library has the flexibility to allow the use of 
pre-existing image data (e.g., phantoms). By means of a series of questions, 
SimSET collects the necessary information to read this data. 

The Photon History Generator is configured using a text file, in which sev¬ 
eral options can be set, like number of decays to simulate, whether to simulate 
SPECT or PET, initial photon energy, simulation time, pointers to the already 
generated activity and attenuation files, effective field of view of the simulation 
etc. 

Once the configuration is done, the simulation can be executed. Several fea¬ 
tures can be retrieved from the photons detected (e.g, binning according to 
energy level, binning according to scattering history, setting of format pre- 
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cision, etc.). In our case, we are interested in the generated sinogram data 
file, which sorts the detected photons based on the transaxial distance and 
azimuthal angle. 

For more details of the SimSET library, the authors suggest the reader the 
citations [16, 49]. 

4.4 Using the NCAT phantom to simulate emission to¬ 
mography in human lungs 

In order to reproduce the breathing movement in emission tomography simu¬ 
lations, a thorax phantom called NCAT phantom has been used. 

NCAT, stands for NURBS-based cardiac torso and was developed by P. Segars 
[41]. It is a model of the human thorax anatomy and physiology created pri¬ 
marily for the nuclear medicine imaging research. The fourth dimensionality 
of the phantom allows modelling of the heart beating and respiratory motions, 
which is of our primary interest. 

Different states within the breathing cycle were acquired and passed as input 
data to the SimSet library. Then, the obtained sinograms were combined to 
simulate movement during the data acquisition. In the next section it is ex¬ 
plained how this respiratory movement information is utilized in the image 
reconstruction algorithm. 
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5 Results 

To test for movement correction, several tests were performed. Starting by a 
simple 2D sphere-shape object, where correction of translation and deforma¬ 
tion of the object was achieved. Followed by 2D and 3D simulations of human 
lungas with simulated respiratory movement. To finalize with a 3D simulation 
of a small lesion added to the lungs to test for movement correction over a 
small structure. 

5.1 Movement correction applied to 2D synthetic images 

Synthetic 2D images simulating a radioactive rod was created, which translates 
a known magnitude. These activity objects were used in a SimSET simulation 
to simulate photons decay. Then, the obtained sinograms (i.e. reference state 
and translated state) were averaged to obtain a final sinogram which simulates 
an instantaneous translation of the radioactive rod. 

Fig. 8 shows the reconstruction obtained by applying MLEM without any 
movement correction. 

From Fig. 8 it is clear to see that MLEM has reconstructed an image com- 



Figure 8: Reconstruction of a moving radioactive rod without movement cor¬ 
rection. 

posed by the two states. Now, as the displacement field is known (a simple 
translation), it is possible to apply the algorithm with movement correction. 
The result is shown in Fig. 9. Fig. 9 shows how the movement correction 
technique has successfully corrected the spatial activity of the simulated rod. 
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Figure 9: Reconstruction of a moving radioactive rod with movement correc¬ 
tion. 



Figure 10: Reconstruction of a linearly deforming radioactive rod without 
movement correction. 

Next, correction of a deforming rod was also achieved. To simulate this situa¬ 
tion, the rod was deformed linearly in one direction, then the same procedure 
applied for the case of the translating rod was applied. Fig. 10 shows the 
reconstruction obtained without correction. 

Fig. 10 shows clearly how the deformation adds to the reconstructed image 
a fuzzy region in the upper border of the rod. Later, the same effect could 
be compared with the fuzzy zone produced in the frontal region of the human 
lungs during simulated breathing. 

The applied displacement vector field can be used to correct for movement, 
whose result is shown in Fig. 11 

The result obtained shows an improvement in the reconstruction, in terms of 
density reduction in pixels where no activity should appear in a free-movement 
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Figure 11: Reconstruction of a moving radioactive rod with movement correc¬ 
tion. 

image reconstruction. Nonetheless, there are some artifacts in the corrected 
image. Our assumption is that during the image reconstruction with move¬ 
ment correction a discretization in the inverse application of the displacement 
vector field caused some pixels to present more density than others. 

To be more realistic, the NCAT phantom was used. First as a 2D phantom 
and then as a 3D one. For the former, two states of the respiratory cycle were 
taken (namely full-inspiration and full-exhalation). Then, SimSET simulations 
were performed and sinograms were combined to simulate breathing during the 
emission tomography exam. For the posterior reconstruction with movement 
correction, the displacement vector field was calculated by using a non-rigid 
registration algorithm called Pasha [6]. 

From Fig. 12 it can be noticed that lungs deformation occurs mainly in the 
upper region, while the deformation in the posterior part is smaller. From the 
combination of both projections sets, simulation of breathing was achieved. 
As in the case of Fig. 10, it is expected that a fuzzy region should appear 
in the frontal region. However, since it is less notorious than the case of 
the deformed radioactive rod, the results in the application of the movement 
correction is better visualized in terms of image-differences between the two 
states of movement used in the simulation. Fig. 13 shows (a)the resulted 
reconstructed image and (b)the result of subtracting the non-corrected image 
and the reference image (i.e., full-expiration). 
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Figure 12: Image reconstruction of (a)full-expiration and (b)full-inspiration 
states. 



(a) Simulation 
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and 




full-expiration. 



Figure 13: (a) Image reconstruction of two combined states of the respiratory 
cycle, (b) Image-difference between non-corrected image and reference image. 
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As it was expected, from Fig. 13, a fuzzy region is visible in the frontal part of 
the lungs. The movement correction was performed by using the displacement 
vector field obtained from the non-rigid registration between the two states of 
movement. Fig. 14 shows the result obtained, the displacement vector field 
used in the correction step and the image difference between the movement 
corrected image and the full-expiration image, which is our reference. From 



(a) Image reconstruc¬ 
tion with movement 
correction of simulated 
breathing. 


(b) Displacement vec¬ 
tor field. 


(c) Image difference of 
corrected image and 
full-expiration state 
(reference). 


Figure 14: (a)Image reconstruction with movement correction of two combined 
states of the respiratory cycle, (b) displacement vector field obtained from 
registration of full-expiration and full-inspiration states, (c)Image difference of 
corrected image and full-expiration state (reference). 


Fig. 14 it can be seen that application of the movement correction technique 
has decreased the fuzzy region produced by the expansion of lungs during 
inspiration. 

To better compare the results, a thresholding can be applied to both image- 
difference images, which is shown in Fig. 15. 
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Figure 15: Comparison of corrected and non-corrected image, by means of 
thresholding image-differences. 
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5.2 Movement correction applied to 3D images 

To test for movement correction in volumetric images, a similar procedure used 
in 2D was performed. First, from the NCAT phantom, a volumetric dataset 
corresponding to the state of full-expiration was extracted. This dataset was 
deformed using a pre-built displacement vector field. This step differs from 
the 2D case, where the displacement vector field used in the correction was 
calculated by using a registration algorithm, which can introduce errors to the 
movement correction step, whereas our primary objective at this point was to 
check the movement correction technique. Fig. 16 shows 3D isosurfaces of full- 
expiration and the deformed volume (rendered with transparency). Besides, 
the displacement vector field has been included as arrow glyphs. 



Figure 16: Full-expiration and deformed volume (in transparent red) isosur¬ 
faces and displacement vector field. 

Each volumetric dataset of dimensions 64x64x64 underwent a SimSET simu¬ 
lation. Then, the obtained sinograms were combined to simulate movement 
during the data acquisition process. 

The MLEM reconstructions were performed with the parallel architecture pre¬ 
sented in section 4.1. In order to compare results in a qualitative way, meshes 
were generated for the reconstructed reference volume (i.e., full-expiration 
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state) and for the non-corrected and corrected reconstructed volumes. Then, 
for each point in the reference surface its corresponding closest point in the 
normal direction, in the non-corrected and corrected surface was found, and 
the distance between these two points was calculated. In this way for the 
non-corrected and corrected images a distance map was obtained. Fig. 17(a) 
shows the distance map obtained between the reference and the non-corrected 
reconstructed image and Fig. 17(b) shows the distance map obtained using 
the movement corrected image. 




(a) non-corrected (b) corrected 

Figure 17: Non-corrected (a) and corrected color distance map (b) obtained 
from calculating the distance between each point in the reference image (full- 
expiration) and the closest point in the normal direction to it. In wireframe 
representation it is shown the non-corrected (a) and corrected (b) meshes. 

From Fig. 17, it can be seen how the distance values between meshes are 
mapped into the original one, going from low values (blue) to high values 
(red). It can be seen from Fig. 17(b) how the movement correction technique 
has produced low distance values (w.r.t the reference image), indicating a bet¬ 
ter image reconstruction. Nonetheless, the corrected image presents a higher 
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distance values at the bottom of the lungs, which is investigated to be a bor¬ 
der condition problem. It is important to remember that since the generated 
meshes come from an iso-surfacing algorithm, the qualitative criteria of the 
correction method will depend on the isovalue selected to create the surfaces. 
However, despite this issue, the qualitative method helps to visualize in some 
degree (depending on the isovalue) the correction achieved by the technique. 
Next, a spherical lesion with a radius of 1cm was added in order to test the 
movement correction over a small structure inside the lungs. Fig. 18 shows 
the activity image for the reference state. 


(a) sagittal plane (b) axial plane (c) coronal plane 

Figure 18: activity image including of lungs and a simulated lesion. The lesion 
is modelled as a sphere with a radius of 1cm. 

The activity images of Fig. 18 were deformed with the displacement vector 
field showed in Fig. 16 to produce a deformed state similar to that of full- 
inspiration. Then, SimSET simulations and combination of sinograms were 
performed to get a sinogram which simulates respiratory movement. This 
sinogram underwent MLEM iterations with and without movement correction. 
To test the results obtained in the movement correction of the lesion, the 
qualitative method, showed previously for the lungs surfaces, was utilized. 
The isosurfaces were generated using a threshold obtained visually from the 
reconstructed images. Fig. 19 shows the distance maps between the non- 
corrected mesh (Fig. 19(a)-(c) wireframe representation) and corrected mesh 
(Fig. 19(b)-(d) wireframe representation) with reference mesh of the lesion. 
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Figure 19: Qualitative test of movement correction for lesion in Fig. 18. 
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As shown in Fig. 19, the movement correction methodology has clearly im¬ 
proved the spatial distribution of the activity in the lesion. However, as it was 
stated, the qualitative method depends on the segmentation of the structure of 
interest. Moreover, as in this case the structure corresponds to a small volume, 
sensitivity in the step of isosurfacing is higher. Therefore, the qualitative test 
is not enough to measure the quality of the results obtained and thus, it should 
be used in conjunction with other method. For this, profiles of intensity were 
generated to compare among reference, non-corrected and corrected profiles. 
Fig. 20 shows profiles of axial slice 43 along its abscise. 

and Fig. 21 shows the evolution of the movement correction along iterations 
of the MLEM algorithm. It can be seen from Fig. 21 how the corrected profile 
gets close to the reference profile with the increase of MLEM iterations. Also, 
it can be noticed that as we get close to the reference (or as the iterations 
proceeds) it is necessary more iterations to improve the correction. This result 
was expected since the technique applies a count-estimate repositioning using 
the movement information, but as the count-estimates are being "improved" 
at each iteration, then the repositioning of these estimations will also improve 
at each iteration. 

Concerning the speed-up obtained by the parallelization of the algorithm and 
the forward and backward accelerators, Fig. 22 shows a plot of speed-up in 
function of the number of processors. 

Figure 23 shows a bar plot of time repartition for the principal inter-iteration 
tasks as communication master-slave, forward projection, backward projection 
and system matrix normalization. For one iteration of a 64 3 volume and a 64 3 
sinogram data. It can be seen from Fig. 23, that the time spend in commu¬ 
nication increase as the number of processors increase, this result is expected 
since a synchronous communication type has been used. Further, the more 
number of slaves are specified, the less time the three main process will take, 
but with the inconvenient that the communication time will increase and it 
will get proportions comparable to those of forward and backward projection 
(e.g,. number of slaves =12 in Fig. 23), which should be readapted if inde¬ 
pendence between communication time and the number of slaves is required. 
However, the selection of the number of slaves should be done considering the 
size of the data to be reconstructed. 
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(a) x-coordinate: 34 


(b) x-coordinate: 35 



(c) x-coordinate: 36 (d) x-coordinate: 37 

Figure 20: Intensity profiles for axial slice 43 around the lesion area. After 
20 MLEM iterations the corrected profiles show a close relationship with the 
reference profiles. 
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Figure 21: Evolution of movement correction in the lesion area. 
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Figure 22: Speed-up in function of the number of processors for the parallel 
implementation of the MLEM with movement correction algorithm 
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Figure 23: Time repartition for communication time, LS process (normaliza¬ 
tion term for each term of the matrix system), FP process (forward projection) 
and BP (backward projection) for a single iteration and in function of the num¬ 
ber of slaves. 
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6 Conclusions 

Induced motion due to patient breathing during an emission tomography exam 
can leads to artifacts in the reconstructed image, artifacts that can produce 
less accurate diagnosis and more important, incorrect radiotherapy planning 
[33] [34] [42], Here, we have presented a methodology to correct for respiratory 
movement in the image reconstruction step. By modifying the way the system 
matrix is calculated, it was possible to us to add information about the res¬ 
piratory movement into the reconstruction algorithm. Furthermore, volume 
deformation of the emissions elements was taken into account, where an emis¬ 
sion element can suffer not only a translation during the breathing cycle but 
also a deformation. 

Considering Figs. 20 and 21, the movement correction is an iterative process, 
which goes attached intrinsically to the evolution of the MLEM algorithm 
itself. In this sense it would be interesting to follow studies in the use of 
weighting schemes to accelerate the movement correction procedure (as it has 
been used for algebraical algorithms like ART and MART). 

The 2-D and 3-D simulations yielded good results, encouraging the application 
of this methodology to real data. However, since the efficacy of the method is 
based on the certitude degree of knowing the movement produced during the 
exam, it will be important to consider a breathing model allowing a good fit 
with the breathing cycle of the patient. In this sense a parametric respiratory 
model is expected to be used. 

Concerning implementational issues, parallelization of the algorithm and ac¬ 
celerators to the steps of forward and backward projection were presented. 
The Bresenham’s algorithm, a classical algorithm in computer graphics, was 
used successfully in the task of accelerating the forward projection step, which, 
until our knowledge, presents a novel contribution to the speed-up of the algo¬ 
rithm. Both, parallelization of the algorithm and the accelerators presented, 
were crucial in the implementation of the MLEM-3D to get iterations times 
of the order of some minutes. This can be an important factor if we think in 
incorporation of the technique in a clinical scenario. 
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